{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import string\n",
    "#https://stackoverflow.com/questions/19726663/how-to-save-the-pandas-dataframe-series-data-as-a-figure\n",
    "import six\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbpath = \"../../fits/\"\n",
    "productpath = \"../../postfit_derivatives/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rois = []\n",
    "df = pd.read_csv(tbpath + 'fit_table_reweighted.csv') \n",
    "rois += list(df.roi.unique())\n",
    "  \n",
    "rois = list(set(rois))\n",
    "rois.remove('US')\n",
    "rois.remove('Gambia')\n",
    "#sort within US and among other coutries then union back\n",
    "roi_us = np.sort([i for i in rois if i[:2]=='US'])[::-1]\n",
    "roi_other = np.sort([i for i in rois if i[:2]!='US'])[::-1]\n",
    "rois = list(roi_us) + list(roi_other)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "[0.00104092 0.01165663 0.02786484 0.05570219 0.15631139]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAK8AAAE+CAYAAAD730AnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARiklEQVR4nO3dfZBVd33H8fcHDEGg1VrWqZABTJxKJq2TNgs1nT63/tNMqtO00ySDD9VkR0dooVqnRqWaWO20QmgwHQs+p4mDKD4wTdvR1FBHk8KmUhNT4gPZmIJGVmNNQhbS8O0f56wcLnf3/u7dXZbv7uc1c+ewv4dzv2fz4fC795zJUURgltGc6S7ArFcOr6Xl8FpaDq+l5fBaWg6vpeXwWlpF4ZX0Zkk7JR2UFJKGenkzSa+Q9BVJT0p6RNL7JfX1si8zlVykkBTAD4D/BC4BfhQRK7p6I2kDsBnYA9wGnAf8GfAQsDoinuiqcpv1SsN7fkQcrP98H7Com/BKWkwV0q8Bl0bE03X75cBngbdExLu6L99ms6Jlw2hwJ+BlwAJg62hw6/3uBg4Caya4f5uFztQHtlX19q42fXcDKyUtOkO12AxxpsK7pN4eatN3CFBjjFmRZ5yh91lQb4+16RtpGXMKSQPAAMDChQsvWbly5eRXZ1PunnvuGY6ISf1m6UyF92i9PRd4sqVvfsuYU0TENmAbQH9/fwwODk5JgTa1JD002fs8U8uGw/V2aZu+pUA0xpgVOVPh3VdvL23T92LggYh4/AzVYjPEpIdX0jJJKyWd02j+DNVyYa2kuY2xlwPnA7dOdh028xWteSW9HFhe/9gHzJP01vrnhyLilsbwjwK/DjwfGAKIiCOS3ga8B/i8pI9RLRfeABwAtkzwOGwWKv3A9hqqQDbdUG/3ALfQQURskvR9YANwE/Aj4OPAX3jJYL0oCm9E/EbpDscbGxEfBj5cui+z8fiWSEvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS6sovJLmSNog6YCkEUkPS9okaWHh/EWSrpN0r6THJA1L+rKkV0nSxA7BZqvSM++NwGbgfmAdsBP4E2C3pHH3Uff/M9WzivdRPSz7ncBc4EPAX/dUuVlEjPsCLgJOAJ9saV8HBHB1h/mX1uNubGmfBxwEftiphtHXJZdcEpYTMBiF/51LXyVn3qsAAVta2rcDR4E1Heb/ZL093PKX5jgwDDxRUIPZaUqe+r6K6sy7t9kYESOS9tf949kL/BB4k6Qh4D+ABcArgUuA13ZZsxlQFt4lwHBEHGvTdwj4ZUnz6jPpaSLiUUm/B7wf+Hij6zHgioj49HhvLmkAGABYtmxZQbk2W5QsGxYA7YILMNIYM57HgfuA9wC/D1wDfBO4TdJLxpsYEdsioj8i+vv6+grKtdmi5Mx7FHjuGH3zG2PakvTzwJeBDRHxvkb7x6gCvV3SBRHxdFnJZpWSM+9hYLGkc9v0LaVaUrRdMtQ2UIV8Z7MxIo4C/wQsB1YUVWvWUBLeffW41c1GSfOBi4HBDvOX1tu5bfqe0bI1K1YS3h1U39Oub2m/lmqte+tog6QLJK1sGXd/vX1Vs1HSs4GXAo9SrX/NutLxjBcR90q6GVgraRdwO3Ah1RW2PcBtjeF3UC0Dmpd8twCvAP66Xv9+CXgOVfifB7ze613rRek/1+uBIaqvrC6juriwFdgYESfGmxgRD0laDWwEfhu4EngS2A+8ISJ29Va6zXZF4a3PjJvq13jjVozR/i2qixJmk8a3RFpaDq+l5fBaWg6vpeXwWloOr6Xl8FpaDq+l5fBaWg6vpeXwWloOr6Xl8FpaDq+l5fBaWg6vpeXwWloOr6Xl8FpaDq+l5fBaWg6vpeXwWloOr6Xl8FpaDq+l5fBaWg6vpeXwWloOr6Xl8FpaDq+l5fBaWg6vpeXwWloOr6Xl8FpaDq+l5fBaWg6vpVUUXklzJG2QdEDSiKSHJW2StLD0jSQ9R9J7JH2z3scRSV+Q9Ku9l2+zWenjW2+ketbwp6iegjn67OFfkPQ7nR7hKmk5cCewCPgA8HXgWcCLOPlUeLOudAyvpIuAdcCuiLii0f4gcBPVs4RvG2P6qH+s3+tFEfGd3ss1O6lk2XAV1VPct7S0bweOAmvGmyzp14BfAf4mIr4j6RxJC3op1qypJLyrgBPA3mZjRIxQPbl9VYf5v1tvvy1pN9UT35+Q9HVJ4wbfbDwl4V0CDEfEsTZ9h4DFkuaNM/+F9XY78Byqp7+/GjgO3CLpj7uo1+zHSj6wLQDaBRdgpDHm+BhjfqLePgb8ZkQcB5D0aeAg8C5JHxnrQ5+kAWAAYNmyZQXl2mxRcuY9Cpw7Rt/8xpixPFlvPzYaXICIeBT4LPAznDw7nyYitkVEf0T09/X1FZRrs0VJeA9TLQ3aBXgp1ZJirLMuwP/U2++26Rv95uGnCuowO0VJePfV41Y3GyXNBy4GBjvMH/2gd16bvtG27xXUYXaKkvDuAAJY39J+LdVa99bRBkkXSFrZMu7TVOvdNZIWNcY+D3gZ8PWI+GYPtdss1/EDW0TcK+lmYK2kXcDtnLzCtodTL1DcASyn+l54dP6jkt4I/ANwt6QPAvOA19XbdZN0LDbLlF4eXg8MUX3qvwwYBrYCGztdGobqQ5ekYeBNwA1U3xvfBVwdEV/qoW4zFBHTXUOx/v7+GBzstMS2s5GkeyKifzL36VsiLS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fSKgqvpDmSNkg6IGlE0sOSNkla2O0bSlog6aCkkPTe7ks2q5SeeW8ENgP3Uz0reCfVs4d3S+r27H090NflHLPTdHz2sKSLqAK7KyKuaLQ/CNwEXMmpD88eb1+/SPUc4zcBm3op2GxUyVnzKqqnuG9pad8OHAXWlLyRpLn1nH8BdnVRo1lbJU99X0X1lPa9zcaIGJG0v+4vsQFYCVzRaaBZiZIz7xJgOCKOtek7BCyWNG+8HUh6PvAO4PqIGOq6SrM2SsK7AGgXXICRxpjxvA84SPWhryuSBiQNSho8cuRIt9NtBisJ71Hg3DH65jfGtCVpDfAS4HUR8VR35UFEbIuI/ojo7+vzlxR2Ukl4D1MtDdoFeCnVkuJ4u4n1nM3A7cB3Jb1A0guA5fWQZ9Vtz+6hdpvlSsK7rx63utkoaT5wMTA4ztxnUn2nexnwjcbrzrp/Tf3zNd0UbQZl3zbsAK6j+n72i432a6nWureONki6ADgnIg7UTU8Af9hmn33A31N9bfYB4KtdV26zXsfwRsS9km4G1kraRbUEuJDqCtseTr1AcQfVkkD13KeAT7TuU9KK+o/fiojT+s1KlJx5oTrrDgEDVEuAYWArsDEiTkxNaWbjKwpvRDxNdTl33Eu6EbGicH9D1Gdns175lkhLy+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkvL4bW0HF5Ly+G1tBxeS8vhtbQcXkurKLyS5kjaIOmApBFJD0vaJGlhwdyflXS9pLslHZH0mKT9kt5SMt9sLKVn3huBzcD9wDpgJ9WDs3dL6rSPVwMbgG8B1wN/DjwAvBP4sqRn9lC3WednD0u6iCqwuyLiikb7g8BNwJWc+uT3Vp8A3h0R/9toe5+kbwBvAV4DvLeH2m2WKznzXkX1kOstLe3bgaPAmvEmR8RgS3BH7ai3P1dQg9lpSsK7CjgB7G02RsQIsL/u78V59faRHufbLFcS3iXAcEQca9N3CFgsaV43byppLvA24P8Yf8mBpAFJg5IGjxw50s3b2AxXEt4FQLvgAow0xnRjC3ApsDEiHhhvYERsi4j+iOjv6+vr8m1sJisJ71Hg3DH65jfGFJF0A7AW2BYR7y6dZ9aqJLyHqZYG7QK8lGpJcbzkzSS9HXgr8CHgtaVFmrVTEt599bjVzUZJ84GLgcGSN6qD+5fAR4BrIiK6qtSsRUl4dwABrG9pv5ZqrXvraIOkCyStbN2BpI1Uwb0FeHVEnOi5YrNax4sUEXGvpJuBtZJ2AbcDF1JdYdvDqd8W3AEsp/peGABJrwfeAXwb+DxwtaTGFB6JiM9N8DhsFuoY3tp6YAgYAC4DhoGtVN8WdDqLjn4PvIxqydBqD+DwWteUaenZ398fg4NFS2w7y0i6JyL6J3OfviXS0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0nJ4LS2H19JyeC2t4vBKmiNpg6QDkkYkPSxpk6SFZ2K+Watuzrw3ApuB+4F1wE6q5w/vllSyn4nONztF0bOHJV1EFbhdEXFFo/1B4CbgSk59gPakzjdrp/SMdxXVk9y3tLRvB44Ca6Z4vtlpSsO7CjgB7G02RsQIsJ+TT3afqvlmpykN7xJgOCKOtek7BCyWNG8K55udpmjNCywA2gUPYKQx5vhkz5c0AAzUPx6TdF/HavNaDAxPdxFT5IWTvcPS8B4FnjtG3/zGmEmfHxHbgG0AkgYjon/8UvOayccnaXCy91m6bDhM9U/7uW36llItCcY6607GfLPTlIZ3Xz12dbNR0nzgYqDT36qJzjc7TWl4dwABrG9pv5ZqrXrraIOkCySt7HV+B9sKx2U1k49v0o9NEVE2UNoKrAU+BdwOXEh1hexLwG9FxIl63BCwPCLUy3yzUt2Edy7VmXMAWEH1qXgHsDEiHm+MG6J9eIvmm5UqDq/Z2WZab4iZyXeqTcKxxRivs+JfKUlvlrRT0sG6rqEe9/MKSV+R9KSkRyS9X1Jf0eSImLYX8HdUH+R2UX142ww8BfwbMGeq55/lxxbAv1Pd99F8/dF0HldLfd8HPgf8ABjqYR8b6v3cSbWcvB54HPgasLDj/Gk8+Iuo7nf4ZEv7uvqArp7K+WfzsTXC8eHpOoaC+s5v/Pm+bsNLdTXxCar7XeY22i+vj/26TvuYzmXDTL5TbdJqkzRP0qJJrG1SRMTBCe7iZVRfk26NiKcb+90NHKTgdzSd4Z3Jd6pNVm1/QBX2xyR9T9JWSc+a1Eqnz+jv4K42fXcDKzv9pZ3O8M7kO9Umo7a9wNupAvxKqrXyWuCLZ+OZuAdL6u2hNn2HqP7lWtKm78dKb8yZCtN2p9oZMOHaIuKXWpo+KumrwF8Bf1pvM1tQb9v9nkZaxrQ1nWfeo0C7G3Wg/E61icyfSlNV299SBf6yXoo6y4wef7vfU9HvaDrDO5PvVJuS2iLiqdF9T7C+s8Hheru0Td9Sqm8cDrfp+7HpDO9MvlNtSmqr558HPDLRAs8C++rtpW36Xgw8EB1uG5jO8J4td6pNhQkdm6SfHmO/N1B9Ttk9eaVOPUnLJK2UdE6j+TPAk8Da+r6X0bGXA+dT8t9vmr/o3srJq1DXAJuorkLdSeMqFDBUldrb/GzHRvX/uLgLeBfwWuCNVN82BNXXSM+czmOra3w58Nb69QjwaOPnl7eMvbOufUVL+xvq9i9QXWF7B9UVtv8GFnWsYZp/AXPrA3iA6lPnIarLqItaxo0V3qL52Y4NeCnwr/WcEaorUfuB64D5031sLYFs97qzJLx136uA/6qP83vAB4HnltTgu8osLf9vliwth9fScngtLYfX0nJ4LS2H19JyeC0th9fScngtLYfX0vp/wrIuwLuynq0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 144x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "font = {'family' : 'sans-serif',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 18}\n",
    "\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "def simpleaxis(ax):\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    \n",
    "\n",
    "def calc_globalquantiles(col):\n",
    "    _025 = 0\n",
    "    _25 = 0\n",
    "    _5 = 0\n",
    "    _75 = 0\n",
    "    _975 = 0\n",
    "    w = 0\n",
    "    for j, roi in enumerate(rois):\n",
    "        try:\n",
    "            sigma2 = (df.loc[(df.roi==roi)&(df['quantile']=='std'), col].values[0])**2\n",
    "            w += 1/sigma2\n",
    "            _025 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.025'), col].values[0]\n",
    "            _25 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.25'), col].values[0]\n",
    "            _5 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.5'), col].values[0]\n",
    "            _75 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.75'), col].values[0]\n",
    "            _975 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.975'), col].values[0]\n",
    "        except:\n",
    "            print()#print(roi)\n",
    "    return _025,_25,_5,_75,_975,w\n",
    "\n",
    "def calc_globalquantiles_searchtheta(theta):\n",
    "    _025 = 0\n",
    "    _25 = 0\n",
    "    _5 = 0\n",
    "    _75 = 0\n",
    "    _975 = 0\n",
    "    w = 0\n",
    "    for j, roi in enumerate(rois):\n",
    "        #find latest week for roi, assume that is post mitigation\n",
    "        for week in range(11,0,-1):\n",
    "            col = theta + \" (week \"+str(week)+\")\"\n",
    "            x = df.loc[(df.roi==roi)&(df['quantile']=='std'), col].values[0]\n",
    "            if np.isfinite(x) and x!='':\n",
    "                break\n",
    "        try:\n",
    "            sigma2 = (df.loc[(df.roi==roi)&(df['quantile']=='std'), col].values[0])**2\n",
    "            w += 1/sigma2\n",
    "            _025 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.025'), col].values[0]\n",
    "            _25 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.25'), col].values[0]\n",
    "            _5 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.5'), col].values[0]\n",
    "            _75 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.75'), col].values[0]\n",
    "            _975 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.975'), col].values[0]\n",
    "#             print(roi)\n",
    "#             print(sigma2)\n",
    "#             print(week)\n",
    "        \n",
    "        except:\n",
    "            print()#print(roi)\n",
    "    return _025,_25,_5,_75,_975,w,week\n",
    "\n",
    "\n",
    "def plot_global(theta,ax,q,pos):\n",
    "    boxes = [\n",
    "            {\n",
    "#             'x': 0,\n",
    "            'fontsize' : 10,\n",
    "            'label' : label_[theta],\n",
    "            'whislo': q[0],    # Bottom whisker position\n",
    "            'q1'    : q[1],   # First quartile (25th percentile)\n",
    "            'med'   : q[2],    # Median         (50th percentile)\n",
    "            'q3'    : q[3],     # Third quartile (75th percentile)\n",
    "            'whishi': q[4],    # Top whisker position\n",
    "            'fliers': []        # Outliers\n",
    "            }\n",
    "        ]\n",
    "    ax.bxp(boxes, widths = 1, \n",
    "            positions=[pos], showfliers=False,  patch_artist=True,#vert=False,\n",
    "               boxprops=dict(facecolor='none',edgecolor='gray'))\n",
    "    return \n",
    "\n",
    "# roi = rois[0]\n",
    "# col = theta\n",
    "# print(df.loc[(df.roi==roi)&(df['quantile']=='0.025'), col])\n",
    "\n",
    "# df['quantile']\n",
    "\n",
    "label_ = {}\n",
    "label_['R0'] = r'R$_{0}$'\n",
    "label_['Rt'] = 'week of April 15th, 2020'\n",
    "label_['car (week 0)'] = 'week 0'\n",
    "label_['car'] = 'week of April 15th, 2020'\n",
    "label_['ifr (week 0)'] = 'week 0'\n",
    "label_['ifr'] = 'week of April 15th, 2020'\n",
    "label_['q'] = ''\n",
    "\n",
    "theta = 'q'\n",
    "# f,ax = plt.subplots(1,1,figsize=(20,5))\n",
    "f,ax = plt.subplots(1,1,figsize=(2,5))\n",
    "\n",
    "_025,_25,_5,_75,_975,w = calc_globalquantiles(theta)\n",
    "q = np.array([_025,_25,_5,_75,_975])/w\n",
    "print(q)\n",
    "# plot_global(theta,ax,q,0)\n",
    "\n",
    "# theta = 'ifr'\n",
    "\n",
    "# _025,_25,_5,_75,_975,w,week = calc_globalquantiles_searchtheta(theta)\n",
    "# q = np.array([_025,_25,_5,_75,_975])/w\n",
    "# print(q)\n",
    "\n",
    "# plot_global(theta,ax,q,1.2)\n",
    "\n",
    "\n",
    "# x = np.array([0.05]+list(np.arange(0.1,1.1,0.1)))\n",
    "# ax.set_xticks(x)\n",
    "# ax2 = ax.twiny()\n",
    "# ax2.set_xticks(x)\n",
    "# ax2.set_xticklabels(np.round(1/x,1))\n",
    "# ax2.set_xlabel('fold difference in unknown vs known cases')\n",
    "# ax.set_xlabel(r'% CAR$_{t}$')\n",
    "# ax.set_xlabel(r'% IFR$_{t}$')\n",
    "# ax.set_xlabel(r'R$_{t}$')\n",
    "# ax.set_xlabel('q')\n",
    "\n",
    "# simpleaxis(ax)\n",
    "# ax.set_xlim((-1,2))\n",
    "# ax.set_ylim((0,0.2))\n",
    "# ax2.set_xlim((0,0.5))\n",
    "# plt.subplots_adjust(left=0.2,top=0.8,bottom=0.2)\n",
    "\n",
    "# plt.subplots_adjust(left=0.4,top=0.9,bottom=0.2)\n",
    "\n",
    "# plt.savefig(productpath + theta + '_barplot.png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
